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! Abstract 

VO . . 

CN) ■ We have investigated the final accretion stage of terrestrial planets from Mars-mass protoplan- 

ets that formed through oligarchic growth in a disk comparable to the minimum mass solar nebula 
(MMSN), through N-body simulation including random torques exerted by disk turbulence due to 
Magneto-Rotat i onal-I nstability. For the torques, we used the semi-analytical formula developed by 
Laughlin et al. ( 2004 ). The damping of orbital eccentricities (in all runs) and type-I migration (in 
some runs) due to the tidal interactions with disk gas are also included. Without any effect of disk 
gas, Earth-mass planets are formed in terrestrial planet regions in a disk comparable to MMSN 
but with too large orbital eccentricities to be consistent with the present eccentricities of Earth 
and Venus in our Solar system. With the eccentricity damping caused by the tidal interaction 
with a remnant gas disk, Earth-mass planets with eccentricities consistent with those of Earth and 
Venus are formed in a limited range of disk gas surface density (~ 10 times MMSN). However, 
in this case, on average, too many (> 6) planets remain in terrestrial planet regions, because the 
damping leads to isolation between the planets. We have carried out a series of N-body simulations 
including the random torques with different disk surface density and strength of turbulence. We 
found that the orbital eccentricities pumped up by the turbulent torques and associated random 
walks in semimajor axes tend to delay isolation of planets, resulting in more coagulation of planets. 
The eccentricities are still damped after planets become isolated. As a result, the number of final 
planets decreases with increase in strength of the turbulence, while Earth-mass planets with small 
eccentricities are still formed. In the case of relatively strong turbulence, the number of final planets 
are 4-5 at 0.5-2AU, which is more consistent with Solar system, for relatively wide range of disk 
surface density (~ 10~ 4 -10 -2 times MMSN). 
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Introduction 



The final stage of terrestrial planet accretion would be coagul ation among protoplanets ( e.g., 



Lissauerlll987l ). The protoplanets form through oligarchic growth (jKokubo k Ida 



1998 



2000) 



so 



that they are called 'oli garchs.' The protop l anets have almost circular orbits initia lly and are iso 
lated from one another (jKokubo k Idalll998l . |2000j) • Their mass is about Mars mass (jKokubo k Ida 



19981 ) in the case of the minimum mass solar nebula (MMSN) model (|Hayashilll98ll ). Long term 



distant perturbations, however, would pump up eccentricities large enough for o rbit crossing, on 



times cales that depend on mass of protoplanets and their orbital separation (jChambers et al 
Because of relatively st rong eccentricity damping due to tidal interaction with a gas disk 



(|Artvmowicall993l ; IWardlll993T). the orbit crossing may not occur until disk gas surface density S g 
decreases below 10 _3 £ gi MMSN (jlwasaki et al.ll2002l ) where S g] MMSN is the surface density of MMSN. 



N-body simulations without any effect of disk gas (e.g. JChambers k Wethe rill 1998; lAgnor k Canup 



19991 ) show that Earth-mass terrestrial planets are formed at ~ 1AU in a disk with a solid sur- 
face density ~ Sd.MMSN as a result of the orbit crossing but with too large orbital eccentricities 
(~ 0.1) to be consisten t with the present eccentricities of Earth and Venus in our Solar system. 
Kominami k Idal (|2002l . 12004 ) performed N-body simulations, taking into account the eccentricity 
damping caused by tidal interaction with a remnant gas disk and found that final eccentricities 
can be small enough to be consistent with those of Earth and Venus. The remnant disk with 



io- 



10 S g ,MMSN allows orbit crossing, but it is still enough to dam p eccentricities of 



Earth-mass planets to < .01 within disk depletion timescale 
20021 : lAenor k Wardl 120021 1. 



10 6 -10 7 years ([Kominami & Ida 



However, iKominami k Idal (|2002l . 12004 ) found that generally > 6 planets remain in terrestrial 
planet region in strong damping cases. If the damping is weaker, number of planets decreases, while 
resultant eccentricities increase. Only in a limited range of the parameters, it sometimes occurs 
that Earth-mass planets with small enough eccentricities (Earth-like planets) are formed and total 
number of formed planets is < 4-5 in a disk similar to MMSN. 



Chambers! (j200ll ) , lO'Brien et al.1 (120061 ) , and iRaymond et al.1 (120061 ) neglected the effects of a 
gas disk in their N-body simulations, but included dynamical friction from remnant planetesimals. 
Although the effect of the dynamical friction is essentially the same as the damping due to a gas 
disk, the number of formed planets is fewer while they have relatively small eccentricities. More 
detailed calculations are needed to clarify the role of dynamical friction from planetesimals. 

Another possibility to reduce the number of formed planets while their eccentricities are kept 
small is a "shaking-up" process to inhibit isolation of planets due to the damping. If gas giants have 
already formed when the orbital crossing starts, e ccentricity excitation by the secular perturbation s 
from the gas giants can provide the shakes (e.g., iLevison k Agnorll2003l ). IKominami k Idal (J2004) 
showed that perturbations from ga s giants in the curren t orbits do not provide enough shakes in 
the presence of disk gas. However, lO'Brien et al.l (|2006l ) reported that the secular perturbations 
from Jupiter and Saturn produce terrestrial planets more consistent with those in our Solar system 
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in the case without the gas damping but with dynamical friction from planetesimals. 



Nagasawa et al.l (|2005l ) considered a passage of a secular resonance during depletion of disk 



gas as a shaking-up mechanism and carried out N-body simulations. v§ resonance passes through 
~ 1AU when S g ~ 10~ 3 -10~ 2 E gi jy[MSN- Hence, after the eccentricity excitation and coagulation of 
protoplanets caused by the resonance passage, the disk gas is still able to damp the eccentricities 
(Eq. [T4"|) . The damping induces inward orbital migration, so that bodies are captured by the 
resonances. The merged terrestrial planets, however, have to be released from the resonances at 
~ 1AU. 

Here, we consider another shaking-up mechani sm, random torques exe r ted by disk turbulence 
due to Magneto-Rotational-Instabili tv fMKD (e.g.. fealbus & HawlevH l99lh. Laughlin etaD (|2004l . 
hereafter referred to as LS A04) and Nelson Sz Papaloizoul ( 2004 ) carried out fluid dynamical simu- 
lations of MRI an d pointed out that the ra ndom torques may significantly influence orbital motions 
of planetesimals. iRice fc Armitagd (|2003l ) studied accretion of a protoplanet taking into account 
the effect of random walks in semimajor axes induced by the random torques and found that the 
random walks expand effective feeding zone of the protoplanet and it may lead to rapid formation 
of a large core of a gas giant. However, since they did not in tegrate o rbits directly, it is not clear if 
their incorporation of random torques is relevant. Actually, iNelsonl (|2005l . hereafter referred to as 
N05) directly integrated orbits of protoplanets in a turbulent disk and found excitation of orbital 
ec centricities as well as ra ndom walks of semimajor axes. The eccentricity excitation was neglected 



Rice fe Armitagd ([2003). Because N05's orbital integration was done simultaneously with fluid 



m 

dynamical simulation of MRI, the orbital integration was limited to 100-150 Keplerian times, which 
is too short to study accretion process of the protoplanets on > 10 6 Keplerian times. 

In order to perform N-body simulations long enough to calculate full stage of accretion of pro- 
toplanets, we adopt the semi-analytical formula for the random torque developed by LSA04 based 
on their fluid dynamical simulations. We directly incorporate the random torques as forces acting 
o n the protoplanets i n the equations of motion. Hence, eccentricity excitation, which was neglected 



m 



Rice fc Armitagd ([2003), is automatically included, as well as a random walk in semimajor axis. 
Because we do not perform fluid dynamical simulation, N-body simulations on timescales ~ 10 7 
years are able to be done. The analytical formula may only roughly mimic the effects of MRI 
turbulence and the method by N05 is more correct. However, our purpose is rather to explore the 
qualitative effects of the turbulence on orbital evolution and accretion of planets on long timescales 
and the dependence on the key parameters on the problem, so that a great quantitative accuracy 
is not important in the present contribution. 

We also include the damping of eccentricities and inclinations direct l y in orbital integrations as 
forces acting on the protoplanets, essentially following Kominami Ida (|2002l . l2004h. However, we 



Kominami et al 



here adopt more exact forms of forces derived bv lTanaka &; Ward! (|2004h (also see 
(|2005l ) and section 2.4). Since in the turbulence, the effect of type-I migration might be greatly 
diminished (N05), we performed both simulations with type-I migration and those without it. When 
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type-I migration is incl uded, we adopt the formula for forces acting on the protoplanets derived by 



Tanaka Ward! (|2004h 



In section 2, we describe the disk model, the formula of forces for the random torques, eccen- 
tricity damping and type-I migration. In section 3.1, we present the results of one planet case in 
order to clearly see the effect of the random torques that we use, on the orbital evolution. The 
results of N-body simulations of accretion of protoplanets that start with 15 protoplanets of 0.2M§ 
are shown in section 3.2. We mainly consider the stage in which disk gas surface density has de- 
clined significantly so that the random torques are relatively weak. However, the weak random 
torques in the stage play important roles to produce terrestrial planets similar to those in our Solar 
system. Section 4 is the conclusion section. 



Model and Calculation Methods 



2.1. Disk model 



Here, we consider a host star with 1M©. Following llda &; Linl (120041 ) . we scale the gas surface 
density S s of disks as 



£g = 2400/ 9 ( I ^r 3 / 2 gcm- 2 , 



(1) 



where f g is a scaling factor; f g = l corresponds to £ g = 1.4S g MMSN- Because current observations 
cannot strictly constrain the radial gradient of E g , we here assume f g is constant with r. Since we 
consider the stage where disk gas has been significantly depleted (f g < 10~ 2 ), optical depth of the 
disk may be low. For si mplicity, we use the temperature distribution in the limit of an optically 



thin disk (H avashilll981 



T = 2.8 x 10 2 (— )~ 1/2 K. 



1AIT 



(2) 



Corresponding sound velocity is 



1-0 x 105(— ) 



-V^cms- 1 . 



(3) 



2.2. Random Torques due to MRI Turbulence 



Turbulence due to Magneto-Rotational Instability (jBalbus Hawlevlll99ll ) is one of candidates 
to account for the observa tionally inferred disk v is cosity, ba s ed on Ha l ine observat i on du e to disk 
accretion onto stars (e.g.. lHartmann et al.lll998l ). iGammid (|1996l ) and ISano et al.l (|2000i ) pointed 
out the existence of "dead zone" around 1AU where the degree of ionization is so small that 
MRI turbulence is suppressed. However, it is not clear that the dead zone exists in the last 
stage of terrestrial planet formation we are considering. If large fraction of dust grains have been 
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transferred into planetesimals in the stage, ionization degree could be high enough that the dead 



zone disappears (jSano et al.ll2000l ). although secondary dust pro duction due to disr uptive collisions 



betwe en planetesimals may a lso be efficient in this stage (e.g., Ilnaba et al.ll2003l ). On the other 



hand. llnutsuka Sand (|2005l ) proposed self-sustained ionization to suggest that dead zone vanishes, 



irrespective of degree of dust depletion. We here assume that disks are MRI-turbulent at ~ 1AU. 

LSA04 modeled the density fluctuations due to the MRI turbulence, based on their fluid 
dynamical simulation. Here we briefly summarize their results with slight modifications. The 
specific force due to density fluctuations exerted on a planet is given by 



tub 



-rv$, (4) 



where 



64£ g r 2 



(5) 

50 



7T 2 M 

$ = 7 r 2 ft 2 ^A c , m , (6) 
i=l 

(r-r e ) 2 _ t 

A c , m = £e ^ cos(mO -<f> c - O c t) sin(7r— ). (7) 

In the above, m is the wavenumber, the dimensionless variable £ has a Gaussian distribution with 
unit width, r and represent the location of the planet in cylindrical coordinates, 0, = y GMq / r 3 
is Keplerian angular velocity at r, r c and <fi c specify the center of the density fluctuation, and O c 
is Q at r c . 7 is the non-dimensional parameter to indicate the strength of the turbulence that we 
introduced instead of LSA04's dimensional parameter A (A = r yr 5 ^ 2 Q 2 ). The pattern speed O c in 
the time-dependent factor allows the mode center to travel along with the Keplerian flow. With 
m specified, the mode extends for a distance 27rr c /m along the azimuthal direction. The radial 
extent is then specified by choosing a = irr c /4:m so that the mode shapes have roughly a 4 : 1 
aspect ratio. The total fluctuation is expressed by superposition of 50 modes at any given time. 
In Eq.©, i in Yli=i expresses individual modes. Each mode comes in and out of existence with 
the time dependence specified above. An individual mode begins at time to an d fades away when 
At > t = t — to- The duration of the mode At is taken to be the sound crossing time of the 
mode along the angular direction, i.e., At = 27rr c /mc s . After the mode has gone, a new mode is 
generated. For the new mode, r c is chosen randomly in the calculation area and <fi is random in 
< 4> < 2-7T. The azimuthal wavenumber m is chosen to be distributed according to a log random 
distribution for wavenumbers in the range 2 < m < 64. 

We modify their formula in two points. The first one is introduction of the non-dimensional 
parameter 7 for strength of turbulence. While LSA04's simulation range was 1.5-3.5 AU, our sim- 
ulations are done around 1AU. Hence, it may be more useful to use the non-dimensional parameter 
than LSA04's dimensional parameter A (A = 7 r 5 / 2 ri 2 ). Since LSA04 used 3.4 AU and 1 year as 
units of length and time and the middle radius of their simulation was 2.5 AU, the values of A in 
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their paper correspond to ~ 1.27(r/lAU) _1// ' 2 . Three dimensional fluid dynamical simulations by 
LSA04 suggest 7 ~ 10 _3 -10 -2 , but the values of 7 may include large uncertainty, so that we explore 
wide range of 7 (also see discussion in section 3.1). The second point is the range of wavenumber 
m. Although 2 < m < 64 in the original formula given in LSA04, inclusion of m = 1 modes could 
be more consistent with global fluid dynamical simulation (G. Laughlin, private communication; 
also see discussion in section 3.1). On the other hand, modes with large m do not contribute to 
orbital changes, because of their short distances (~ a = irr c /4m) for effective forces and rapid time 
variations with timescale 2Trr c /mc s . Hence, in order to save calculation time, we cut off high m 
modes with m > 6, that is, A cm is set to be zero when m > 6, in the summation of 50 modes in 
Eq. ([6]) The torque exerted on the planet with mass M is 

1 (9$ 

rtub = r x -— x r x M (8) 
r ou 

50 

= - 7 rMr 2 n 2 ^mA Sim , (9) 
i=i 



where 



(r-r p f . , n , „ ~ x „ t 



A s , m = £e ° A sin(m6> - <f> c - Q c i) sin(7r— ). (10) 

Figure [l^, shows the scaled random torques, J2i=i m ^-s,m, with 2 < m < 64. In Fig. [lb, m > 6 
modes are cutted off. Low frequency patterns, which contribute to orbital evolution, are similar 
between these results. Actually, we found through orbital calculations like in Fig. [2] that m < 5 are 
enough to reproduce orbital evolution with fully counting all m modes. If m = 1 mode is included, 
the amplitude of the random torques does not change, but low frequency patterns change. We will 
present the results of N-body simulations with 2 < m < 5 in section 3.2, but we also carried out 
calculations with 1 < m < 5 and will discuss the effects of m = 1 modes. 

[Figure [Q 



2.3. Secular Torques due to Disk-Planet Interactions 



As shown below, the random torques given by Eq. ([9]) induce random walks of semimajor axes 
of planets and pump up their orbital eccentricities. Tidal interactions with a laminar disk monoton- 
ically decreases the semimajor axes and damp the eccent ricities (and the inclinations). The secula r 



inward migration is known as "type- 1 migration" (e.g., IWard 



1986 



19971 : iTanaka & Wardl 12002 1 . 



Since mean flow in turbulent disks coincides with flow in laminar disks, interactions with the mean 
flow may induce the secular orbital migration and eccentricity damping even in turbulent disks. In 
turbulent disks, however, N05 reported that type-I migration might be greatly diminished while 
the eccentricity damping still works. Non-linear effects associated with the random fluctuations 
(e.g., the temporary activation of corotation torques or temporary disruption of the pressure buffer) 
could be responsible for the slowing down. Alternatively, relatively high eccentricities excited by 
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the random torques, which is > h/r ~ 0.05 obtained by N05 where /i is disk scale height, could 
affect the type-I migration (e.g.. iPapaloizou &: Larwwodll2000l ). In our N-body simulations, eccen- 
tricities are also pumped up to > 0.05 by perturbations among protoplanets except for the last 
phases well after orbital crossing. Hence, we performed two series of simulations: one is without 
type-I migration and the other is with it. 

We summarize the secular changes in laminar disks below. Both torques from inner and outer 
disks damp orbital eccentricities and inclinations, since the gravitational interactions with disk gas 
causes similar effect of dynamical friction. The damping timescales are (jTanaka &; Wardll2004l ) 



^damp,e 


e 
e 


^damp 

0.78 


^damp,i 


i 


^damp 

0.54 


^damp 




£ g a 2 




= 240/; 





V K 
a ,2 



1AU' 



years, 



(11) 
(12) 

(13) 

(14) 



where a is the semimajor axis of the planet and t>K is the Keplerian velocity at a. 



On the other hand, the torque from an inner disk increases semimajor axis, while tha t from 
an outer disk decreases it. Since the outer torque is generally greater than the inner one ([Wardl 
198a : iTanaka Wardl 120021 ) , the torque imbalance induces inward migration (type-I migration) . 
For the radial gradient of S g cx a -1,5 , the torque imbalance, which is negative definite, is given by 



(ITanaka Wardl 12002^ 



r mig = -2.17 A 2 (^) 2 £ g a^ 2 . 

IVIq C s 



Migration timescale due to this torque is 

(l/2)Mfia 2 



"mig 



a 
a 



i mig 

0.23(— r x (- 



Mr. 







5.0 xl0 4 (— )- x ( 

Ma) 



)-l(±L)2 n -l 
VK 

Y^f^fg 1 years. 



(15) 

(16) 
(17) 
(18) 



2.4. Orbital Integration 



We integrate orbits of 15 protoplanets with O.2M0 that initially have orb it s of s mall e and % 
(~ 0.01) with separation 6?~h, following initial conditions in iKominami fc Idal (|2002), where Hill 
radius rjj is defined by 



(— ) 



1/3, 



0.007( 



M 
0.2M f 



1/3, 



(19) 
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Initial angular distributions are set to be random. Calcu lation starts from the phase when the 
orbital crossing starts. The result of iKokubo k, Idal (|2000l ) shows that the eccentricities of proto- 
planets produced through oligarchic growth are about ~ 10 -3 , so that the protoplanets are well 
isolated. However, the protoplanets will eventually start orbital crossing by long- term mutual dis 



tant perturbations on a ti me scale depending on their orbital separation, mass (jChambers et al 



1996 ), initial eccent ricities (jYoshinaga et al.lll999i ). and how much gas is around the protoplanets 
(jlwasaki et al. 2002 ). Since we are concerned with orbital crossing stage, we start the calculation 
with relatively high eccentricities e = 10~ 2 , supposing the eccentricities have already increased and 
orbital crossing is ready to start. The initial inclinations are also set to be i = 10~ 2 . 

The basic equations of motion of particle k at in heliocentric coordinates are 



d 2 r k 
dt 2 



-GM C . 



\ r k\' c 



"I'-^damp ~H -f tub ~)~ -^migi 



rk 



r k 



3 



3 I on . |3 



(20) 



where k, j = 1, 2, 15, the first term is gravitational force of the central star, the second term is 
mutual gravity between the bodies, and the third term is the indirect term. -Fdamp and -F m i g are 
specific forces for the damping of eccentricities and inclinations and type-I migration, and -Ftub is 
specific force due to the turbulence (Eq. 0]). Their detailed expressions are described below. Note 
that in our simulations, mass of bodies is larger than 0.2M m , so that aerodynamical drag forces are 
neglected compared with Fd amp and -F m i g (e.g.. IWardlll993l ). 



We integrate orbits with the fourth-order Hermite scheme. When protoplanets collide, perfect 
accretion is assumed. After the collision, a new body is created, conserving total mass and mo- 
mentum of the two colliding protoplanets. The physical radius of a protoplanet is determined by 
its mass and internal density as 

/ 3 M M/3 
47T pp 

The internal density pp is set to be 3 gem -3 . 



(21) 



Tanaka &: Wardl (|2004l ) derived, through three-dimensional linear analysis, 



^damp,r 



F d . 



amp, z 



mig,r 



mig,i 

F ■ 

1 mig,. 



0<-> 4 <¥ 

M Q c s M q 



Mq C s Mq 



)n(2A c r [v e - rn\ + A s r v r ) 
rO] + A a g v r ) 



( 



M 



)(^) 4 ( 



)n(A c z v z + A s z zU) 



Mq C 



'' K ) 2 ( 



Mq 



)9,v K 



0, 



(22) 

(23) 

(24) 
(25) 
(26) 
(27) 
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where 



At = 0.057 Af, = 0.176 
-0.868 A s e = 0.325 
-1.088 At = -0.871. 



A1 



Note that there is a typos in Fd am p,z in iTanaka Ward (120041 ) . The factor (2A c z v z + A s z zQ) 
should be {A'jv? + A% z tt) as in Eq. (|24p . Note also that the other factors in the expressions in 
Kominami et al.l ([20051 ) have minor typos; the above expressions are correct ones. Eccentricities 
are damped by Fd amPjr and Fd amPi 0, while inclinations are damped by -Fdamp^- Semimajor axes 
are decreased by F mig Q (= T m ; g /Mr) where a and r are identified because of small e and i. The 
evolution of e, i and a by orbital integration of one body with the above forces completely agrees 
with the analytically derived evolution with Eqs. (|14p and (|18p . 



The force due to turbulence, F tu b 



-rV<3? (Eq. H]), is given by 

50 



-^tub,r 

Ftuhfi 
-Ptub,2 



7 rr^ 2 ( 1 + 
i=l ^ 



2r(r 



(7 



50 



jTrQ 2 mA s 
i=i 

0, 



where A c m and A StTn are defined by Eqs. ([7]) and ([10 



(28) 

(29) 
(30) 



As seen above, -Fdamp and F m \g are parameterized by the disk surface scaling factor fg for 
given M and r of planets, and -F^ub by f g and 7 (r oc f g ). Therefore, f g and 7 are parameters for 
our calculations. As discussed in section 1, we will consider the stages in which disk gas has been 
significantly depleted, so that the cases of f g = 10 -4 and 10 -2 are mainly studied. Although the 
most likely value of 7 might be ~ 10 _3 -10 -2 , it would include large uncertainty (also see discussion 
in section 3.1), so that the cases of 7 = lO^lO" 1 andl are studied. For comparison, non-turbulent 
(7 = 0) cases are also calculated. 



3. Results 

3.1. One Planet Case 

To see the effects of the turbulent forces given by Eqs. (|28|) and ([29]) on orbital changes and how 
they depend on f g and 7, we first carry out simulations with one planet embedded in a turbulent 
disk. Figures [2] show evolution of semimajor axis a and orbital eccentricity e of a planet of O.2M0 
obtained by orbital integration with .Ftub i n the case of 7 = 10 _1 and f g = 10~ 2 . The initial 
conditions are a = 1AU and e = 0. -Fdamp and -^mig are not included. As expected, a random walk 
of a and excitation of e are observed. 
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[FigureE] 
[Figure [3] 

In order to quantify the random walks, we performed 100 similar runs with different random 
numbers for the random torques, but still using 7 = 10 _1 and f g = 10~ 2 . At each time, the distri- 
butions of deviation in semimajor axis Aa from the initial position (1AU) and orbital eccentricity e 
for the 100 runs are fitted as Gaussian distributions to obtain the standard deviations as functions 
of time. Hereafter, the standard deviations are also denoted by Aa and e. Figures [3] show the 
evolution of Aa and e obtained by the numerical calculations. The evolution curves are fitted as 

Aa - 1.8 x 1(T 6 (— ^— ) 1 / 2 AU, (31) 
lyear 

e ~ 2.7 x lQ-Sf-J— ) 1 ' 2 . (32) 
lyear 

near 1AU. The dependence of i 1 / 2 would reflect diffusion characteristics. (If -Fdamp is included, e 
approaches an equilibrium value.) We have carried out the same procedures for 7 = 10~ 2 , 10 _1 and 
f g = 10~ 2 , 10 _1 to derive the dependence of 7 and f g as 

Aa ~ 2 x l(T 3 /o7(— ^— ) 1/2 AU. (33) 
lyear 

e ~ 3 x 10- 2 / 5 7(t^) 1/2 . (34) 
lyear 

Note that the random walks are independent of planet mass M. In the above calculations, we used 
m = 2-5 modes. With m = 2-64, we obtained very simular results. 

Equations (|33p and (|34p give Aa and e that are 10-100 times smaller than those obtained 
by LSA04 and N05. We found that the analytically modeled random torques almost cancel out 
in time and the net change is only ~ 0.001 of total change for m = 2-5. Since both LSA04 and 
N05 used global fluid codes to follow orbits of protoplanets, m = 1 modes might be included. 
Since m = 1 modes have the longest duration and distance for effective force, it would induce 
asymmetry between positive and negative torques to produce larger Aa and e. We have carried 
out similar calculations, including m = 1 modes and found that Aa and e are 10 times larger than 
Eqs. (|33p and (|34p . The inclusion of m = 1 modes may mostly resolve the difference from the 
results by LSA04 and N05, but the approximated semi-analytical torque formula could be still too 
symmetric, compared with the global fluid dynamical simulations. Hence, the results of 7 = lO^ 1 
and 1 (with m = 2-5 modes), which are larger than the numerically inferred value 7 ~ 10" 3 -10- 2 , 
are also pertinent for the evolution of planets in realistic turbulent disks. In the N-body simulations 
shown in section 3.2, perturbations from other protoplanets also induce some asymmetry and Aa 
and e may be much larger than Eqs. (133p and (|34p during the period in which protoplanets undergo 
relatively close encounters. 
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If type-I migration works on the timescale given by Eq. <\18\i. the migration length near 1AU 

is 

A^xlO-V^X^AU. (35) 



From Eqs. f|33|) and (|35p . it is expected that if 

^ 3xl0V( 0^ r2yearS ' (36) 

type-I migration will dominate over the random walk. Figured] shows the evolution of the semimajor 
axis with both effects of type-I migration and turbulent fluctuations of 7 = 0.1. The planet starts 
secular inward migration after t ~ 3 x 10 3 years, consistent with the above estimate. Note, however, 
that in the turbulent disks, it is not clear that type-I migration speed is still the same as that 
predicted by the linear calculation (N05). 

[Figure H] 



3.2. Accretion of Protoplanets in a Turbulent Disk 



Because we will compare the results with lKominami &: Ida! (|2002l ) and because type-I migration 



might be greatly diminished in turbulent disks, in many runs we calculate accretion and the orbital 
evolution of protoplanets in a turbulent disk without the effect of type-I migration. We carry out 
simulations with various f g and 7. We denote a run with f g = 10"°, 7 = 10 -/3 as RUNo^, where 
k (k = a, b, c) represent different initial angular distribution of the protoplanets. In some runs, the 
effect of type-I migration is included, which we denote as RUN a^ a f. Table 1 shows simulation 
parameters for individual runs with 7 < 1 and m = 2-5 (28 runs). Two runs were carried out with 
7 = 10 (m = 2-5). We also carried out 18 runs with inclusion of m = 1 modes and found that 
slightly smaller 7 produce similar results to the cases without m = 1 modes. To avoid confusion, 
we will only present the detailed results with m = 2-5. 



3.2.1. Case with f g = 10~ 2 

First we show the results with f g = 10~ 2 . The orbital evolution of RUN2oo a , RUN2 3a , RUN2i a , 
and RUN2oa are shown in Figs. [5h., b, c, and d, respectively. The thick solid lines represent 
semimajor axes a. The thin dashed lines represent pericenters a(l — e) and apocenters a(l + e). 
Thicker solid lines represents more massive planets. With f g = 10 -2 , the damping time scale 
Tdamp ^ 1-2 x 10 5 years for M = 0.2M e . 



[Figure [5] 
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Since RUN2oo a does not include the random t orques (7 = 0), the evolution in Fig. [5^ is very 



similar to that shown by iKominami Idal (|2002l ). In this case, a planet of O.6M0 with small 



eccentricity (~ 0.0001) is formed. However, global orbital crossing lasts for only ~ 5 x 10 5 years 
because of the rather strong eccentricity damping. Consequently, the number of surviving planets 
are 8, which is much greater than that in the present Solar system. The runs with very weak 
turbulence of 7 = 10 -3 in Fig. [5b shows a similar result to the non-turbulent case. For 7 = and 
10 -3 , the number of surviving planets is always 8 or 9 (Table 1). 

The effects of turbulence are pronounced in the cases of 7 = 10 _1 (Fig. [Sb) and 7 = 1 (Fig. EH). 
The random walk and eccentricity excitation induced by the turbulence tend to inhibit isolation 
of the planets. In Fig. \5jp (RUN2i a ), the duration of orbit crossing is elongated, while 8 planets 
still survive. (The same number of planets survive also in RUN2k, and RUN2i c .) According to 
the eccentricity excitation effect, the eccentricities of final planets are slightly larger than in the 
previous two cases, however, they are still smaller than the present free eccentricities of Earth and 
Venus, because the damping that increases with planet mass eventually overwhelms the turbulent 
excitation that is independent of the planet mass. In Fig. [5H (RUN2o a with 7 = 1), the large 
random walk enhances the number of collision events (10 events), so that the number of surviving 
planets drastically decreases to 4. In RUN2 f, and RUN2o c , the number of surviving planets is also 
4 or 5. 

In Fig.[5H (7 = 1), secular inward migration is found, although type-I migration is not included. 
This migration is induced by the damping of eccentricities that are continuously excited by the 
random torques, since orbital angular momentum, ^ GM®a{l — e 2 ), is almost conserved during 
the eccentricities damping. In the run with extremely large 7 (= 10), the turbulent excitation is so 
strong that all the planets are removed from terrestrial planet region by the inward migration. 



3.2.2. Case with f g 



10' 



The evolution in severely depleted disks with f g = 10 , RUN4oo a , RUN4 3a , RUN4i a , and 
RUN4oa, are shown in Figs.[6K, b, c, and d, respectively. f g = 10~ 4 corresponds to Td amp ~ 1.2 x 10 7 
years for M = 0.2M^. RUN4 ooa shows the result without the effect of turbulence. The weak 
damping due to the small surface density of a gas disk elongates the period during which the 
eccentricities are high enough to allow orbital crossing (~ 1 x 10 7 years). As a result, a larger 
planet (M = 1.4M®) than in RUN2oo a is formed. Since e is damped down to ~ 0.01, this planet 
is very similar to Earth. However, the number of surviving p lanets is 6 (Fig. EK), wh ich is larger 
than that in the present Solar system, as is the case shown bv lKominami Idal (|2002l ). The mean 
number of surviving planets of RUN4oo a , RUN4oob and RUN4oo C is 6.3. 



[Figure [6] 
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In the turbulence cases, the mean number of surviving planets is 5.8 (7 = 10~ 3 ), 5.7 (7 = 1CT 1 ), 
and 4.7 (7 = 1). As the turbulence becomes stronger, the number of final planets decreases. In 
RUN4i a with 7 = lO^ 1 , global orbital crossing lasts on more than 10 7 years (Fig. [6k), while 
RUN43a does not show such clear elongation of orbital crossing (Fig. [6)3). The turbulent excitation 
for eccentricities is still weaker than the tidal damping for Earth-mass planets as long as 7 < 1, so 
that their final eccentricities are still < 0.01. (In the run with extremely large 7 (= 10), we found 
that the eccentricities are not sufficiently damped.) In general, probability for final planets to be 
similar to present terrestrial planets in our Solar system is larger for f g = 10 -4 than for f g = 10 -2 . 

The accretion timescales in weak turbulence cases a re a few times 10 6 ye ars a fter orbit crossing 
starts . Those in strong turbulence cases are ~ 10 7 years. Ilwasaki et al.l (|2002l ) and lKominami Idal 
((2002) suggested that orbit crossing does not start until f g decays down to ~ 10~ 3 . If the effect 



of turbulence is taken into account, orbit crossing may start at the stage of larger f g . If the 
condition of f g < 10~ 3 is applied and exponential decay from initial f g ~ 1 with decay timescale 
Tdep is assumed, the orbit crossing starts at ~ 7rd ep ~ 10 7 -10 8 years. Thus, the total accretion 
timescales in the present model ar e not in contradiction to the Earth formation age inferred from 
Hf- W chronology ~ 4 x 10 7 years 



Yin and Ozima 


2003; 


Kleine et al. 


2002, 


2004) 



If disk depletion is only due to viscous diffusion, it may be possible that such small-mass 
remnant disks remain for 10 7 -10 8 years. However, if disk dispersal due to stellar EUV is efficient, 
it would be difficult to preserve such small-mass remnant disks. Spitzer surve y found that 25% o f 
B-A members and 10% of F-K members in Pleiades cluster show IR excess (jNadya et al.ll2006l ). 
The excess might imply the existence of small-mass remnant gas disks, but it might also be due 
to secondary dust generation in gas free environments. Observation of gas components for clusters 
at 10 7 -10 8 years is needed to examine the role of the tidal damping due to remnant disks on final 
orbital configuration of terrestrial planets. 



3.2.3. Eccentricities and Feeding Zones 

Once orbit crossing starts, the velocity dispersion is pumped up to surface escape velocity v esc 
of planets by close encounters. Corresponding eccentricity is given by 

e ~ V — = 0.34(-^)V6 ( _^)V3 ( « (37) 
vk 3gcm 6 M® 1AU 

Figure [7] shows eccentricity evolution of all bodies in RUN2oo a (non-turbulent case) and RUN2o a 
(strongly turbulent case with 7 = 1). During orbit crossing (t < 1 x 10 6 years), the mean ec- 
centricities are almost same in the two cases. Since eccentricity excitation due to random torques 
evaluated by Eq. (|34p is significantly smaller than Eq. (|37p . the eccentricities during orbit crossing 
are mostly determined by mutual planetary perturbations. 



[Figure [7] 
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In RUN2oo a , global orbit crossing ceases at t > 1 x 10 6 years and then the eccentricities are 
secularly decreased by the damping due to -Fdamp- However, in RUN2o a , close encounters still 
occasionally occur at t > 1 x 10 6 years, so that the eccentricity damping is slower. Even with 
relatively strong turbulence of 7 = 1 of this run, the tidal damping of eccentricities overcomes the 
turbulent excitation for Earth-mass planets. (But, this is not the case for 7 = 10.) 

For 7 = 1 and f g = 10 -2 , diffusion length due to random torques is evaluated by Eq. ([33]) as 
Aa ~ 10 _2 (i/10 6 year) 1 / 2 AU, which is much smaller than orbital separation among the protoplan- 
ets. However, as suggested before, planetary perturbations may inhibit cancellation of the torques 
and induce much larger Aa (and e). Furthermore, even if the turbulence itself does not directly 
expand the feeding zones of the planets, scattering by close encounters among protoplanets that 
are induced by the random torques allows the feeding zones to effectively expand. In the N-body 
simulations, the two effects are indistinguishable. 



3.2.4- Effects of Type I Migration 



Here, the results with type-I migration (calculations with -F m ig) are shown, although it is not 
clear that type-I migration actually operates in turbulent disks (N05). In RUTN^ooa/, f g = 10 -2 and 
the turbulenc e is not includ e d (7 = 0). More systemat i c inve stigations in the non-turbulence cases 
were done by iMcNeil et all (|2005l ) and baisaka et al.1 (fcood ). Although [McNeil et all (j2005T ) and 
Daisaka et al.l (|2006l ) included the effects of small planetesimals as well, RUN2 ooa / shows similar 
properties to their calculations (Fig. [8^): planets in inner regions tend to fall onto the host stars 
while those in outer regions could survive. Since collision events are limited by the loss of inner 
planets, the mass of the largest surviving planets is 0.4M®. (Inclusion of protoplanets in more 
outer region might increase the mass of final planets.) 



[Figure [8] 



Figure[8b and[8b show RUN2i a / of 7 = 10 _1 and RUN2o a / of 7 = 1. Even with relatively strong 
turbulence, the tendency to migrate inward does not change, compared with the non-turbulent 
case in Fig. [Hh.. Equation (|36p shows that type-I migration is dominant over the random walk after 
t ~ 3 x 10 5 7 2 (M/0.2M®) -2 years, so that the random walk cannot halt the inward migration on 
timescales ~ 10 6 years. The inward migration is rather accelerated by the damping of eccentricities 
that are continuously excited by the random torques (section 3.2.1). 

For larger f g , since the random torques are stronger, the accelerated migration is more pro- 
nounced. Thus, our results suggest that random migration superposed to type-I migration would 
not be able to solve the problem that planets tend to be lost from the terrestrial planet region. 
The problem can be solved only if the turbulent fluctuations somehow inhibit the underlying type-I 
migration, as found by N05, or if planets at ~ 1AU are formed by surviving protoplanets originally 
at > 1 AU. 
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4. Conclusions 



We have investigated the final accretion stage of terrestrial planets from Mars-mass protoplan- 
ets in turbulent disks, through N-body simulation. Gravitational interactions with gas disks exert 
the following three effects on the protoplanet orbits: 

1. damping of eccentricities e and inclinations i, 

2. type-I migration (secular decrease of semimajor axis a), 

3. random-walks of a and stochastic excitation of e. 



The effect 3) has not been included in N- body simulations o f plane t accretion in the previous works. 
We adopt the same simulation setting of [Kominami &: Idal (|2002j ) that included only the effect 1): 
initially 15 protoplanets of 0.2M^ are set with orbital separations of several Hill radii in terrestrial 
planet regions, corresponding to MMSN. In our N-body simulations, the effects 1) and 3) were 
included. The effect 2) was examined in section 3.2.4. We incorporated random torques exerted 
by disk turbulence due to MRI as forces directly acting on protoplanets in the equations of motion 
for orbital integration. We adopted the semi-analytical formula for t he random torque s devel oped 
by LSA04 with slight modifications. Compared with the results of IKominami 8z Idal (|2002l ). we 
investigated the effects of disk turbulence on planet accretion. 

The past N-body simulations neglecting the gas disk showed that the coagulation between 
protoplanets result in planets of about Earth-mass but with the eccentricities higher than the 
present terrestrial planets in our Solar system. If the effect 1) is included, when £ g ~ 10~ 4 - 
10 _3 E gi MMSN, the damping allows initiation of orbit crossing to form an Earth-mass planet(s), 
while it damps the eccentricities sufficiently after p lanets are isolated. Ho wever, > 6 planets tend 
to remain, because of isolation due to the damping (|Kominami Idall2002l ). (Note that in the case 
of damping by dyna mical friction from remnant planetesimals the number of planets is reduced 
(jO'Brien et al.ll2006l )). We found that the newly incorporated effect 3) tends to inhibit isolation of 
planets, resulting in more coagulations of planets, while the eccentricity damping is still effective. 
As a result, 4-5 planets with small eccentricities are formed in relatively wide parameter range: 
gas surface density S g ~ 10~ 4 -10 _2 £ gi MMSN ; and MRI turbulence strength 7 ~ 10" 1 -! (slightly 
smaller 7 if m = 1 modes of density fluctuation are included). 

LSA04's prescription for the random torques that we adopted has highly symmetric properties 
and the exerted torques almost completely cancel out in time averaging. As a result, the diffusion 
length and eccentricity excitation obtained by one planet calculations are generally too small to 
play a direct role in expanding feeding zones of protoplanets in late phase in which disk gas is 
significantly depleted (Eqs. [331 and I34p . However, planetary perturbations may break the symme- 
try. Furthermore, in more realistic turbulence, the torques may include m = 1 modes and be less 
symmetric. These effects inhibit the torque cancellation to induce much larger Aa and e. Further- 
more, even if the enhanced effects are still too small to directly expand the feeding zones, such 
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small effects can be enough to break the isolation of the protoplanets, thus allowing them to have 
distant encounters with each other. The encounters in turn induce larger random oscillations of 
the semimajor axes, effectively enhancing the feeding zone of each planet. We also found through 
calculations with all the effects 1), 2) and 3) that the random walks do not decelerate (rather 
accelerate) the type-I migration, although it is not clear that type-I migration actually operates in 
turbulent disks. 

Although the prescription for the random torques would include large uncertainty, we have 
demonstrated that the random torques tend to decrease number of final planets while they keep 
formation of Earth-mass planets with small eccentricities, which is more consistent with the present 
Solar system. Since the random torques are independent of mass of bodies, small planetesimals also 
suffer the random torques. N-body simulations starting from smaller planetesimals in turbulent 
disks will be presented in a separate paper. 
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Table 1: Initial parameters and final results for each run. M\ is the mass of the largest planet in 
final state. e\ is the time averaged eccentricity of the largest planet, taken after its isolation takes 
place. 
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Fig. 1. — An example of the random torques that we adopted. The (scaled) total torque mK s 
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Fig. 4. — Orbital evolution with both effects of type-I migration and turbulent fluctuations with 
f g = 1(T 2 and 7 = 10 _1 . 
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Thp results nf N-hnrlv simiilatirms with f_ = 1 fT 4 fat -v = fRTTN4„J CM -v = in -3 
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Thp results nf N-hnHv simnla.t.inns with f' = 
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